﻿/*
    Copyright 2010, 2011 Eric Wong
	contact@optimalcycling.com
	http://www.optimalcycling.com
	Optimal Cycling - Advanced power pacing program for cyclists
	
	This file is part of Optimal Cycling.

    Optimal Cycling is free software: you can redistribute it and/or modify
    it under the terms of the GNU Lesser General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Optimal Cycling is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public License
    along with Optimal Cycling.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef __OPENCL_VERSION__
	#include "Kernels.h"
#endif

#ifdef _USRDLL
//extern "C"
//{
	__declspec(dllexport)
#endif
__kernel void
OCSimRun(
	const int currSolutionMode, 
    const int whichArray,
	const sizeT numIndividuals,
	const sizeT numSimPoints,
	  
	__global real* powerValA, 
	__global real* powerValB,
	__global const real2* powerLimits,
	const int powerAllSame, 

	__global real* fitness,
	__global sizeT* correspondingIndex,
	__global sizeT* bestFitIndex,
	const sizeT numIndividPerPop, 
	const sizeT numPopulations,

	__global const real4* coursePosPoints,
	__global const real* coursePosLinear,
	const real coursePosLinearMax,
	__global const real4* coursePosDirection,
	__global const real4* coursePosInfo,
	__global const real4* riderProp,
	__global const real4* wind, 
	__global const real4* climate,

	__global const real* dragCoeffWheelYawAngles,
	__global const real2* dragCoeffWheels, 
	const sizeT dragCoeffWheelsNumPts,

	__global const ModelInfo* modelInfo, 
	__global const real4* iConstants,

	__global uint4* randNumGenSeeds,
	__global const SolverInfo* solverInfo,
	__global SolverVars* solverVars,
	__global const PowerMetricsInfo* powerMetricsInfo,
	__global const real2* powerMetricsCriticals,

	__global real* tempWorkingSpace,

	__global real* bestTime,
	__global real8* bestPowerResults,
	__global real8* bestKinematics,

	__global real2* fitStats)
{
	__global real* powerVal;
	__global real* powerValSecond;

	#ifdef _OPENMP
		int nIndex;
	#else
		sizeT nIndex;
	#endif	

	if (whichArray == 0)
	{
		powerVal = powerValA;
		powerValSecond = powerValB;
	}
	else
	{
		powerVal = powerValB;
		powerValSecond = powerValA;
	}

	switch(currSolutionMode)
	{
		case SolutionMode_Initialize:
			if (get_global_id(0) == 0)
			{			
				fitStats[0].s0 = LARGE_REAL;
				fitStats[0].s1 = 0;

				bestTime[0] = LARGE_REAL;
				solverVars->currGen = 0;
			}

			#ifdef _OPENMP
				omp_set_num_threads(OMP_NUM_THREADS);
				#pragma omp parallel for
			#endif
				for (nIndex = get_global_id(0); nIndex < numIndividuals; nIndex += get_global_size(0)) 
				{ 
					sizeT upperPowerIndex = correspondingIndex[nIndex] + numSimPoints;
					sizeT counter = 0;
					sizeT individIndexInPop;
					sizeT i;

					fitness[nIndex] = LARGE_REAL;
										
					// Set heuristic seeds and initialize all others with random values
					individIndexInPop = nIndex % numIndividPerPop;
					if (powerAllSame != 0)	// Genes must be all same in an individual
					{
						if (individIndexInPop < solverInfo->heuristicSeedIndexLow)	// Best guess provided by user
						{							
							real steadyPower = clamp(powerValA[correspondingIndex[nIndex]], powerLimits[0].s0, powerLimits[0].s1);
							for(i = correspondingIndex[nIndex]; i < upperPowerIndex; i++)
							{
								powerValA[i] = steadyPower;
								powerValB[i] = powerValA[i];								
							}
						}
						else	// Individuals who are randomly seeded; Heuristic seeding not done, to be implemented?
						{
							real steadyPower = randRealBounded(&randNumGenSeeds[nIndex], powerLimits[0].s0, powerLimits[0].s1);
							for(i = correspondingIndex[nIndex]; i < upperPowerIndex; i++)
							{
								powerValA[i] = steadyPower;
								powerValB[i] = powerValA[i];
							}
						}
					}
					else
					{						
						if (individIndexInPop < solverInfo->heuristicSeedIndexLow)	// Best guess provided by user
						{
							for(i = correspondingIndex[nIndex]; i < upperPowerIndex; i++)
							{
								powerValA[i] = clamp(powerValA[i], powerLimits[counter].s0, powerLimits[counter].s1);
								powerValB[i] = powerValA[i];
								counter++;
							}
						}
						else if (individIndexInPop >= solverInfo->heuristicSeedIndexLow && individIndexInPop <= solverInfo->heuristicSeedIndexHigh)	// Individuals we heuristically seed
						{
							real steadyPower = randRealBounded(&randNumGenSeeds[nIndex], (real)100, (real)350);
							for(i = correspondingIndex[nIndex]; i < upperPowerIndex; i++)
							{
								powerValA[i] = clamp(steadyPower, powerLimits[counter].s0, powerLimits[counter].s1);
								powerValB[i] = powerValA[i];
								counter++;
							}
						}
						else	// Individuals who are randomly seeded
						{
							for(i = correspondingIndex[nIndex]; i < upperPowerIndex; i++)
							{
								powerValA[i] = randRealBounded(&randNumGenSeeds[nIndex], powerLimits[counter].s0, powerLimits[counter].s1);
								powerValB[i] = powerValA[i];
								counter++;
							}
						}
					}
				}
			break;


		// Randomly move individuals from one population to another
		case SolutionMode_Migrate:
			if (get_global_id(0) == 0)			
			{
				sizeT individ1, individ2, i;

				if (numPopulations <= 1)
					return;

				for (i = 0; i < solverInfo->numIndividMigrate; i++)
				{
					individ1 = randULongBounded(&randNumGenSeeds[0], 0, numIndividuals - 1);
					individ2 = randULongBounded(&randNumGenSeeds[0], 0, numIndividuals - 1);

					if (individ1 != individ2)
						SwapEx(fitness, correspondingIndex, individ1, individ2);
				}
			}
			break;

		case SolutionMode_OutputAll:
			if (get_global_id(0) == 0)
			{	
                __global real* currPowerVals = &powerVal[correspondingIndex[bestFitIndex[0]]];
				bestTime[0] = ObjectiveFunc(&tempWorkingSpace[0], currPowerVals, powerLimits, iConstants, solverVars->maxTime,
						solverVars, powerMetricsInfo, powerMetricsCriticals, KERNEL_VARS_TO_MODEL, bestPowerResults, bestKinematics);
			}
			break;


		case SolutionMode_EvaluateFitness:
			#ifdef _OPENMP
				omp_set_num_threads(OMP_NUM_THREADS);
				#pragma omp parallel for
			#endif
				for (nIndex = get_global_id(0); nIndex < numIndividuals; nIndex += get_global_size(0)) 
				{ 
					__global real* currPowerVals;

					currPowerVals = &powerVal[correspondingIndex[nIndex]];

					fitness[nIndex] = ObjectiveFunc(&tempWorkingSpace[nIndex * numSimPoints], currPowerVals, powerLimits, 
						iConstants, solverVars->maxTime, solverVars, powerMetricsInfo, powerMetricsCriticals, KERNEL_VARS_TO_MODEL, 0, 0);
				}
			break;

		// Using quick sort to sort by increasing fitness for each population
		// Only works for positive values.
		case SolutionMode_SortByFitness:
			#ifdef _OPENMP
				omp_set_num_threads(OMP_NUM_THREADS);
				#pragma omp parallel for
			#endif
				for (nIndex = get_global_id(0); nIndex < numPopulations; nIndex += get_global_size(0)) 
				{
					sizeT popStartIndex = nIndex * numIndividPerPop;
					QuickSortEx(&fitness[popStartIndex], &correspondingIndex[popStartIndex], numIndividPerPop);
				}
			break;
		
		case SolutionMode_NoOptimize:
			if (get_global_id(0) == 0)
			{	
                __global real* currPowerVals = &powerVal[0];

				fitStats[0].s0 = LARGE_REAL;
				fitStats[0].s1 = 0;

				solverVars->currGen = 0;
				fitness[0] = ObjectiveFunc(&tempWorkingSpace[0], currPowerVals, powerLimits, iConstants, solverVars->maxTime,
						solverVars, powerMetricsInfo, powerMetricsCriticals, KERNEL_VARS_TO_MODEL, bestPowerResults, bestKinematics);
				bestTime[0] = fitness[0];
			}
			break;

		// Updates the generation count and finds the fitness individual in the populations
		// Must have sorted populations before calling
		case SolutionMode_IncrementParams:
			if (get_global_id(0) == 0)
			{								
				real bestFitVal = fitness[0];
				real fitValTemp;
				sizeT bestFitIndexTemp;
				sizeT fitCount = 0;
				bestFitIndex[0] = 0;

				for (nIndex = 1; nIndex < numPopulations; nIndex++) 
				{
					bestFitIndexTemp = nIndex * numIndividPerPop;
					fitValTemp = fitness[bestFitIndexTemp];
					if (fitValTemp < bestFitVal)
					{
						bestFitVal = fitValTemp;
						bestFitIndex[0] = bestFitIndexTemp;
					}
				}

				bestTime[0] = bestFitVal;

				fitValTemp = 0;				
				for (nIndex = 0; nIndex < numIndividuals; nIndex++) 
				{
					if (fitness[nIndex] < solverVars->maxTime)
						fitCount++;
					fitValTemp += fitness[nIndex];
				}
				fitStats[0].s0 = fitValTemp / numIndividuals;
				fitStats[0].s1 = fitCount;

				solverVars->currGen = solverVars->currGen + 1;
			}
			break;

		// This is the actual solver logic that searches out the optimal solution
		case SolutionMode_DiffEvo:
			#ifdef _OPENMP
				omp_set_num_threads(OMP_NUM_THREADS);
				#pragma omp parallel for
			#endif
				for (nIndex = get_global_id(0); nIndex < numIndividuals; nIndex += get_global_size(0)) 
				{ 	
					real testcost, F, Cr, Sr, currMaxTime;
					sizeT rnd0, rnd1;
					sizeT currCorrespIndex = correspondingIndex[nIndex];
					__global real* childPowerGenes = &powerValSecond[currCorrespIndex];
					__global real* parentPowerGenes = &powerVal[currCorrespIndex];
					__global real* powerTempSpace = &tempWorkingSpace[nIndex * numSimPoints];

					sizeT popNumber = nIndex / numIndividPerPop;
					//sizeT individIndexInPop = nIndex % numIndividPerPop;

					__global real* bestGenesThisPop = &powerVal[correspondingIndex[popNumber * numIndividPerPop]];
					
					sizeT counter = 0;
					sizeT stratIndexer = randULongBounded(&randNumGenSeeds[nIndex], 0, numSimPoints - 1);												

					F = randRealBounded(&randNumGenSeeds[nIndex], solverInfo->diffEvoFLimits.s0, solverInfo->diffEvoFLimits.s1);
					Cr = randRealBounded(&randNumGenSeeds[nIndex], solverInfo->diffEvoCrLimits.s0, solverInfo->diffEvoCrLimits.s1);							

					ArrayCopy(parentPowerGenes, childPowerGenes, numSimPoints);

					do 
					{
						rnd0 = randULongBounded(&randNumGenSeeds[nIndex], 0, numIndividPerPop - 1);
					} while (rnd0 == nIndex);

					do 
					{
						rnd1 = randULongBounded(&randNumGenSeeds[nIndex], 0, numIndividPerPop - 1);
					} while ((rnd1 == nIndex) || (rnd1 == rnd0));

					//DERandToBest1Exp
					do
					{
						// The jittering improves convergence. Keep small.
						childPowerGenes[stratIndexer] += (F * (1.0f + 0.001f*(randRealTraditional(&randNumGenSeeds[nIndex]) - 0.5f)))
						  * ((bestGenesThisPop[stratIndexer] - parentPowerGenes[stratIndexer]) 
						  + (powerVal[correspondingIndex[rnd0] + stratIndexer] - powerVal[correspondingIndex[rnd1] + stratIndexer]));			

						childPowerGenes[stratIndexer] = clamp(childPowerGenes[stratIndexer], powerLimits[stratIndexer].s0, powerLimits[stratIndexer].s1);						  

						stratIndexer = ++stratIndexer % numSimPoints;
					} while ((randRealTraditional(&randNumGenSeeds[nIndex]) < Cr) && (++counter < numSimPoints));

					
					// Flat power smoothing
					if (numSimPoints >= 20 && randRealTraditional(&randNumGenSeeds[nIndex]) < 0.1f)
					{
						sizeT smoothNum = randULongBounded(&randNumGenSeeds[nIndex], 2, numSimPoints / 10);						
						real smoothingSum = 0;
						real firstSmoothVal;
						int validFlatPower = 1;
						stratIndexer = randULongBounded(&randNumGenSeeds[nIndex], 0, numSimPoints - smoothNum);
						
						firstSmoothVal = childPowerGenes[stratIndexer];
						for (counter = 0; counter < smoothNum; counter++)
						{			
							if (fabs(childPowerGenes[stratIndexer + counter] - firstSmoothVal) > 20)
							{
								validFlatPower = 0;
								break;
							}
							smoothingSum += childPowerGenes[stratIndexer + counter];
						}

						if (validFlatPower)
						{
							smoothingSum = smoothingSum / smoothNum;

							for (counter = 0; counter < smoothNum; counter++)
							{			
								childPowerGenes[stratIndexer + counter] = clamp(smoothingSum, 
									powerLimits[stratIndexer + counter].s0, powerLimits[stratIndexer + counter].s1);
							}
						}
					}

					// 3-Point Smoothing
					if (randRealTraditional(&randNumGenSeeds[nIndex]) < 0.5f)
					{
						Sr = randRealBounded(&randNumGenSeeds[nIndex], solverInfo->diffEvoSrLimits.s0, solverInfo->diffEvoSrLimits.s1);
						for (counter = randULongBounded(&randNumGenSeeds[nIndex], 1, numSimPoints - 1); counter < numSimPoints - 1; counter++)
						{			
							if (randRealTraditional(&randNumGenSeeds[nIndex]) < Sr)
							{
								childPowerGenes[counter] = clamp((childPowerGenes[counter - 1] + childPowerGenes[counter] + childPowerGenes[counter + 1]) / 3, 
									powerLimits[counter].s0, powerLimits[counter].s1);
							}
							else
								break;
						}
					}

					// Butterworth smoothing
					if (randRealTraditional(&randNumGenSeeds[nIndex]) < 0.5f)
					{
						real bwCutOffF = randRealBounded(&randNumGenSeeds[nIndex], solverInfo->bwCutOffF.s0, solverInfo->bwCutOffF.s1);
						real bwDataCollectF = ONE;
						sizeT smoothNum = randULongBounded(&randNumGenSeeds[nIndex], 0, numSimPoints);

						if (smoothNum >= 5)
						{
							sizeT startSmoothIndex = randULongBounded(&randNumGenSeeds[nIndex], 0, numSimPoints - smoothNum);
							BWSmoothData(bwCutOffF, bwDataCollectF, smoothNum, startSmoothIndex, powerTempSpace, childPowerGenes, powerLimits, 1);
						}
					}

					if (solverVars->raceTimeMeasureStartIndex == 0)
						currMaxTime = min(solverVars->maxTime, fitness[nIndex]);
					else
						currMaxTime = solverVars->maxTime;

					testcost = ObjectiveFunc(&tempWorkingSpace[nIndex * numSimPoints], childPowerGenes, powerLimits, 
						iConstants, currMaxTime, solverVars, powerMetricsInfo, powerMetricsCriticals, KERNEL_VARS_TO_MODEL, 0, 0);

					if (testcost <= fitness[nIndex])
					{						
						fitness[nIndex] = testcost;	
					}
					else	
					{
						ArrayCopy(parentPowerGenes, childPowerGenes, numSimPoints); // Propagate the old
					}

				}
			break;
	}
}

#ifdef _USRDLL
//}
#endif





